/*
 *  GaussSeidel.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 3/1/11.
 *  Copyright 2011 ICS - USI - Lugano. All rights reserved.
 *
 */

#ifndef NSE_SOLVER_GaussSeidel_H
#define NSE_SOLVER_GaussSeidel_H


#include "Solver.h"
/**
 * \brief The GaussSeidel class ... TODO Doc
 */ 

namespace nsesolver {


class GaussSeidel : public Solver< GaussSeidel >{
public:
	
	/**
	 * The empty Constructor
	 */
	GaussSeidel() {}
	
	template< class MatrixType >
	void step(const MatrixType &A, MatrixType &u, const MatrixType &f) const
	{
		typedef typename MatrixType::Scalar Scalar;
		
		for (int i = 0; i < u.size(); i++) {
			Scalar sum = 0;
			for (int j = 0; j < i; j++) {
				sum += A(i,j) * u[j];
			}
			
			for (int j = i + 1; j < u.size(); j++) {
				sum += A(i,j) * u[j];
			}
			
			u[i] = -(sum - f[i])/A(i,i);
		}
	}
	
	
	template< class MatrixType >
	void stepSparse(const MatrixType &A, MatrixType &u, const MatrixType &f, const CRS< Array< int > > &crs) const
	{
		typedef typename MatrixType::Scalar Scalar;
		
		int accumlatedJ = 0;
		for (int i = 0; i < u.size(); i++) {
			Scalar sum = 0;			
			for (int j = 0; j < crs.getRowCount(i); j++, accumlatedJ++) {
				if (crs.getColumnIndex(accumlatedJ) == i) {
					continue;
				}
				sum += A(i, crs.getColumnIndex(accumlatedJ)) * u[crs.getColumnIndex(accumlatedJ)];
			}
			u[i] = -(sum - f[i])/A(i,i);
		}
	}
	
	
	/**
	 * The Destructor
	 */
	~GaussSeidel()
	{}
	
};
	
}

#endif
